Clustering Samples with vst

K-means Clustering

No recognisable best k, so chose k=5

set.seed(123)
wss = sapply(1:15, function(k) kmeans(datExpr_redDim, k, nstart=25)$tot.withinss)
plot(wss, type='b', main='K-Means Clustering')
best_k = 3
abline(v = best_k, col='blue')

datExpr_k_means = kmeans(datExpr_redDim, best_k, nstart=25)
clusterings[['km']] = datExpr_k_means$cluster

Hierarchical Clustering

Chose k=7 as best number of clusters.

Clusters seem to be able to separate ASD and control samples pretty well and there are no noticeable patterns regarding sex, age or brain region in any cluster.

Younger ASD samples seem to be more similar to control samples than older ASD samples (pink cluster has most of the youngest samples), perhaps oder samples were only diagnosed in extreme cases and now milder ASD cases get diagnosed a well and milder cases have milder gene expression differences? The yellow cluster is made of young ASD samples. Most old ASD samples are also close together.

Colors:

• Diagnosis: Blue=control, Green=ASD

• Sex: Pink=Female, Blue=Male

• Brain region: Pink=Frontal, Green=Temporal, Blue=Parietal, Purple=Occipital

• Age: Purple=youngest, Yellow=oldest

h_clusts = datExpr_redDim %>% dist %>% hclust %>% as.dendrogram
# h_clusts %>% plot
best_k = 7
clusterings[['hc']] = cutree(h_clusts, best_k)

create_viridis_dict = function(){
  min_age = datMeta_redDim$Age %>% min
  max_age = datMeta_redDim$Age %>% max
  viridis_age_cols = viridis(max_age - min_age + 1)
  names(viridis_age_cols) = seq(min_age, max_age)
  
  return(viridis_age_cols)
}
viridis_age_cols = create_viridis_dict()

dend_meta = datMeta_redDim[match(gsub('X','',labels(h_clusts)), rownames(datMeta_redDim)),] %>% 
            mutate('Diagnosis' = ifelse(Diagnosis_=='CTL','#008080','#86b300'), # Blue control, Green ASD
                   'Sex' = ifelse(Sex=='F','#ff6666','#008ae6'),                # Pink Female, Blue Male
                   'Region' = case_when(Brain_lobe=='Frontal'~'#F8766D',        # ggplot defaults for 4 colours
                                        Brain_lobe=='Temporal'~'#7CAE00',
                                        Brain_lobe=='Parietal'~'#00BFC4',
                                        Brain_lobe=='Occipital'~'#C77CFF'),
                   'Age' = viridis_age_cols[as.character(Age)]) %>%            # Purple: young, Yellow: old
            dplyr::select(Age, Region, Sex, Diagnosis)
h_clusts %>% set('labels', rep('', nrow(datMeta_redDim))) %>% set('branches_k_color', k=best_k) %>% plot
colored_bars(colors=dend_meta)

Consensus Clustering

Samples are grouped into two big clusters with 8 and 7 subclusters, respectively.

*Output plots in clustering_samples_04_05 folder

Independent Component Analysis

Following this paper’s guidelines:

0. Run PCA and keep enough components to explain 60% of the variance

1. Run ICA with that same number of nbComp as principal components kept to then filter them

2. Select components with kurtosis > 3

3. Assign obs to genes with FDR<0.01 using the fdrtool package

## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr
## 
## Step 1... determine cutoff point
## Step 2... estimate parameters of null distribution and eta0
## Step 3... compute p-values and estimate empirical PDF/CDF
## Step 4... compute q-values and local fdr

Not a good method for clustering samples because:

1. ICA does not perform well with small samples (see Figure 4 of this paper)

2. Warnings: (Warning in fdrtool(x, plot = F): There may be too few input test statistics for reliable FDR calculations!)

3. Leaves almost half of the observations (59) without a cluster:

ICA_clusters %>% rowSums %>% table

## .
##  0  1  2  3 
## 59 24  2  1

ICA_clusters %>% mutate(cl_sum=rowSums(.)) %>% as.matrix %>% melt %>% ggplot(aes(Var2,Var1)) + 
  geom_tile(aes(fill=value)) + xlab('Clusters') + ylab('Samples') + 
  theme_minimal() + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + coord_flip()

WGCNA

SFT.R.sq starts in 0.83 and then 0.846 but it then decreases and never achieves a higher value, so chose power=2

best_power = datExpr_redDim %>% t %>% pickSoftThreshold(powerVector = seq(1, 30, by=2))

##    Power SFT.R.sq   slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.8380  1.6600          0.812  43.400    45.600  58.30
## 2      3   0.0135  0.0343         -0.254  19.500    20.100  35.80
## 3      5   0.2460 -0.3180          0.149  11.300    10.700  25.30
## 4      7   0.6310 -0.5260          0.624   7.320     6.310  19.10
## 5      9   0.6590 -0.6920          0.616   5.110     3.930  15.10
## 6     11   0.6950 -0.7960          0.735   3.760     2.640  12.30
## 7     13   0.7230 -0.7810          0.761   2.870     1.770  10.20
## 8     15   0.8140 -0.8490          0.886   2.260     1.260   8.59
## 9     17   0.8430 -0.8920          0.885   1.820     0.962   7.35
## 10    19   0.7800 -0.9520          0.782   1.490     0.690   6.36
## 11    21   0.8040 -0.9870          0.799   1.240     0.510   5.55
## 12    23   0.7840 -1.0100          0.772   1.050     0.404   4.88
## 13    25   0.7840 -0.9680          0.774   0.898     0.333   4.31
## 14    27   0.8300 -0.9280          0.871   0.775     0.278   3.84
## 15    29   0.8090 -0.9530          0.811   0.676     0.234   3.43

network = datExpr_redDim %>% t %>% blockwiseModules(power=2, numericLabels=TRUE)
clusterings[['WGCNA']] = network$colors
names(clusterings[['WGCNA']]) = rownames(datExpr_redDim)

It finds two main clusters and leaves only 3 observations without a cluster

table(clusterings[['WGCNA']])

## 
##  0  1  2 
##  3 50 33

Gaussian Mixture Models with hard thresholding

Points don’t seem to follow a Gaussian distribution no matter the number of clusters, chose 4 points following the best k from K-means because the methods are similar

n_clust = datExpr_redDim %>% Optimal_Clusters_GMM(max_clusters=80, criterion='BIC', plot_data=FALSE)
plot(n_clust, type='l', main='Bayesian Information Criterion to choose number of clusters')

best_k = 3 # copying k-means best_k
gmm = datExpr_redDim %>% GMM(best_k)
clusterings[['GMM']] = gmm$Log_likelihood %>% apply(1, which.max)

Plot of clusters with their centroids in gray

gmm_points = rbind(datExpr_redDim, setNames(data.frame(gmm$centroids), names(datExpr_redDim)))
gmm_labels = c(clusterings[['GMM']], rep(NA, best_k)) %>% as.factor
ggplotly(gmm_points %>% ggplot(aes(x=PC1, y=PC2, color=gmm_labels)) + geom_point() + theme_minimal())

rm(wss, datExpr_k_means, h_clusts, cc_output, cc_output_c1, cc_output_c2, best_k, ICA_output, 
   ICA_clusters_names, signals_w_kurtosis, n_clust, gmm, gmm_points, gmm_labels, network, dend_meta, 
   best_power, c, viridis_age_cols, create_viridis_dict)